home *** CD-ROM | disk | FTP | other *** search
/ Super PC 34 / Super PC 34 (Shareware).iso / spc / UTIL / DJGPP2 / V2 / DJTST200.ZIP / tests / libc / ansi / math / cephes / mtst.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-13  |  11.7 KB  |  497 lines
  1. /*   mtst.c
  2.  Consistency tests for math functions.
  3.  Inverse function tests and Wronksians for random arguments.
  4.  With NTRIALS=10000, the following are typical results for
  5.  IEEE double precision arithmetic:
  6.  
  7. Consistency test of math functions.
  8. Max and rms relative errors for 10000 random arguments.
  9. x =   cbrt(   cube(x) ):  max = 1.54E-016   rms = 6.50E-018
  10. x =   atan(    tan(x) ):  max = 2.83E-016   rms = 7.68E-017
  11. x =    sin(   asin(x) ):  max = 4.44E-016   rms = 7.01E-017
  12. x =   sqrt( square(x) ):  max = 1.57E-016   rms = 2.70E-017
  13. x =    log(    exp(x) ):  max = 1.88E-014   rms = 2.02E-016
  14. x =   tanh(  atanh(x) ):  max = 4.00E-016   rms = 8.29E-017
  15. x =  asinh(   sinh(x) ):  max = 2.04E-016   rms = 1.56E-017
  16. x =  acosh(   cosh(x) ):  max = 2.67E-014   rms = 2.95E-016
  17. x = pow( pow(x,a),1/a ):  max = 3.38E-012   rms = 3.39E-014
  18. Legendre  ellpk,  ellpe:  max = 1.55E-011   rms = 1.74E-013
  19. x =  ndtri(   ndtr(x) ):  max = 4.09E-014   rms = 6.40E-016
  20. Absolute error criterion (but relative if >1):
  21. lgam(x) = log(gamma(x)):  max = 4.49E-016   rms = 8.80E-017
  22. x =    cos(   acos(x) ):  max = 4.44E-016   rms = 1.00E-016
  23. Absolute error and only 2000 trials:
  24. Wronksian of   Yn,   Jn:  max = 7.12E-011   rms = 1.59E-012
  25. Wronksian of   Kn,   Iv:  max = 7.07E-012   rms = 5.25E-013
  26.  
  27. */
  28.  
  29. /*
  30. Cephes Math Library Release 2.1:  December, 1988
  31. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  32. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  33. f*/
  34.  
  35. #include <stdlib.h>
  36. #include <stdio.h>
  37. #include <math.h>
  38. #include <string.h>
  39. #include <signal.h>
  40. #ifdef  LDOUBLE
  41. #include <mathext.h>
  42. #endif
  43. #include "mconf.h"
  44.  
  45. #define NTRIALS 10000
  46. #define WTRIALS (NTRIALS/5)
  47. #define STRTST 0
  48.  
  49. #ifdef LDOUBLE
  50. #define double long double
  51. #undef fabs
  52. #define fabs fabsl
  53. #define sqrt sqrtl
  54. #define cbrt cbrtl
  55. #define exp expl
  56. #define log logl
  57. #define tan tanl
  58. #define atan atanl
  59. #define sin sinl
  60. #define asin asinl
  61. #define cos cosl
  62. #define acos acosl
  63. #define pow powl
  64. #define tanh tanhl
  65. #define atanh atanhl
  66. #define sinh sinhl
  67. #define asinh asinhl
  68. #define cosh coshl
  69. #define acosh acoshl
  70. #define drand drandl
  71. #define pow2 pow2l
  72. #define log2 log2l
  73. #define pow10 pow10l
  74. #define log10 log10l
  75. #define MINLOG MINLOGL
  76. #define MAXLOG MAXLOGL
  77. #define LOG2E  LOG2EL
  78. #define CPI    PIL
  79. #define PIO2   PIO2L
  80. #define GFORMAT "%Lg"
  81. #define P18EFORMAT "%.20LE"
  82. #define P4EFORMAT "%.4LE"
  83. #define P2EFORMAT "%.2LE"
  84. #else
  85. #define GFORMAT "%g"
  86. #define P18EFORMAT "%.18E"
  87. #define P4EFORMAT "%.4E"
  88. #define P2EFORMAT "%.2E"
  89. #ifdef  __DJGPP__
  90. #include <float.h>
  91. #define _set_cw87(x)    _control87(x, 0xffffffffU)
  92. #endif
  93. #endif
  94.  
  95. /*
  96. double          ndtr(), ndtri(), ellpe(), ellpk(), gamma();
  97. double          fabs(), sqrt();
  98. double          cbrt(), exp(), log(), tan(), atan();
  99. double          sin(), asin(), cos(), acos(), pow();
  100. double          tanh(), atanh(), sinh(), asinh(), cosh(), acosh();
  101. double          jn(), yn(), iv(), kn();
  102. */
  103.  
  104. /* Provide inverses for square root and cube root: */
  105. double
  106. square(x)
  107.     double          x;
  108. {
  109.     return (x * x);
  110. }
  111.  
  112. double
  113. cube(x)
  114.     double          x;
  115. {
  116.     return (x * x * x);
  117. }
  118.  
  119. /* lookup table for each function */
  120. struct fundef
  121. {
  122.     char           *nam1;    /* the function */
  123.     double          (*name) ();
  124.     char           *nam2;    /* its inverse  */
  125.     double          (*inv) ();
  126.     int             nargs;    /* number of function arguments */
  127.     int             tstyp;    /* type code of the function */
  128.     long            ctrl;    /* relative error flag */
  129.     double          arg1w;    /* width of domain for 1st arg */
  130.     double          arg1l;    /* lower bound domain 1st arg */
  131.     long            arg1f;    /* flags, e.g. integer arg */
  132.     double          arg2w;    /* same info for args 2, 3, 4 */
  133.     double          arg2l;
  134.     long            arg2f;
  135. /*
  136.     double arg3w;
  137.     double arg3l;
  138.     long arg3f;
  139.     double arg4w;
  140.     double arg4l;
  141.     long arg4f;
  142. */
  143. };
  144.  
  145.  
  146. /* fundef.ctrl bits: */
  147. #define RELERR 1
  148. #define EXPSCAL 4
  149.  
  150. /* fundef.tstyp  test types: */
  151. #define POWER 1
  152. #define ELLIP 2
  153. #define GAMMA 3
  154. #define WRONK1 4
  155. #define WRONK2 5
  156. #define WRONK3 6
  157.  
  158. /* fundef.argNf  argument flag bits: */
  159. #define INT 2
  160.  
  161. extern double   MINLOG;
  162. extern double   MAXLOG;
  163. extern double   LOG2E;
  164. extern double   CPI;
  165. extern double   PIO2;
  166.  
  167. /*
  168. define MINLOG -170.0
  169. define MAXLOG +170.0
  170. define PI 3.14159265358979323846
  171. define PIO2 1.570796326794896619
  172. */
  173.  
  174. /*
  175. #ifndef LDOUBLE
  176. #define NTESTS 13
  177. #else
  178. #define NTESTS 11
  179. #endif
  180. */
  181.  
  182. struct fundef   defs[] =
  183. {
  184.     {"   tan", tan, "  atan", atan, 1, 0, 1, 0.0, 0.0, 0,
  185.      0.0, 0.0, 0},
  186.     {"  asin", asin, "   sin", sin, 1, 0, 1, 2.0, -1.0, 0,
  187.      0.0, 0.0, 0},
  188.     {"square", square, "  sqrt", sqrt, 1, 0, 1, 170.0, -85.0, EXPSCAL,
  189.      0.0, 0.0, 0},
  190.     {"   exp", exp, "   log", log, 1, 0, 1, 340.0, -170.0, 0,
  191.      0.0, 0.0, 0},
  192.     {" atanh", atanh, "  tanh", tanh, 1, 0, 1, 2.0, -1.0, 0,
  193.      0.0, 0.0, 0},
  194.     {"  sinh", sinh, " asinh", asinh, 1, 0, 1, 340.0, 0.0, 0,
  195.      0.0, 0.0, 0},
  196.     {"  cosh", cosh, " acosh", acosh, 1, 0, 1, 340.0, 0.0, 0,
  197.      0.0, 0.0, 0},
  198.     {"pow", pow, "pow", pow, 2, POWER, 1, 25.0, 0.0, 0,
  199.      50.0, -25.0, 0},
  200.     {"  exp2", pow2, "  log2", log2, 1, 0, 1, 340.0, -170.0, 0,
  201.      0.0, 0.0, 0},
  202.     {" exp10", pow10, " log10", log10, 1, 0, 1, 340.0, -170.0, 0,
  203.      0.0, 0.0, 0},
  204. #if 0
  205.     {" ellpe",  ellpe,   " ellpk",  ellpk, 1, ELLIP, 1, 1.0,    0.0,  0,
  206.      0.0, 0.0, 0},
  207.     {"  ndtr",   ndtr,   " ndtri",  ndtri, 1, 0, 1,   10.0,   -10.0,   0,
  208.      0.0, 0.0, 0},
  209. #endif
  210.     {"  acos", acos, "   cos", cos, 1, 0, 0, 2.0, -1.0, 0,
  211.      0.0, 0.0, 0},
  212. #if 0
  213.     {"  Iv",     iv,   "  Kn",     kn, 2, WRONK2, 0,  30.0, 0.1,  0,
  214.      20.0, 0.0, INT},
  215. #endif
  216.     {NULL, NULL, NULL, NULL, 0, 0, 0,0.0, 0.0, 0, 0.0, 0.0, 0}
  217. };
  218.  
  219. static char    *headrs[] =
  220. {
  221.     "x = %s( %s(x) ): ",
  222.     "x = %s( %s(x,a),1/a ): ",    /* power */
  223.     "Legendre %s, %s: ",    /* ellip */
  224.     "%s(x) = log(%s(x)): ",    /* gamma */
  225.     "Wronksian of %s, %s: ",
  226.     "Wronksian of %s, %s: ",
  227.     "Wronksian of %s, %s: "
  228. };
  229.  
  230. static double   cy1 = 0.0;
  231. static double   y2 = 0.0;
  232. static double   y3 = 0.0;
  233. static double   y4 = 0.0;
  234. static double   a = 0.0;
  235. static double   x = 0.0;
  236. static double   y = 0.0;
  237. static double   z = 0.0;
  238. static double   e = 0.0;
  239. static double   max = 0.0;
  240. static double   rmsa = 0.0;
  241. static double   rms = 0.0;
  242. static double   ave = 0.0;
  243.  
  244. void mysig(int sig);
  245. int
  246. main()
  247. {
  248.     double          (*fun) ();
  249.     double          (*ifun) ();
  250.     struct fundef  *d;
  251.     int             i, k=0, itst;
  252.     int             m, ntr;
  253.     _set_cw87(0x137f);
  254.     ntr = NTRIALS;
  255.     printf("Consistency test of math functions.\n");
  256.     printf("Max and rms relative errors for %d random arguments.\n",
  257.        ntr);
  258.  
  259. /* Initialize machine dependent parameters: */
  260. #ifdef __TURBOC__
  261. #define FACT 0.99
  262. #else
  263. #define FACT 1.0000
  264. #endif
  265.     defs[0].arg1w = CPI;                /* tan */
  266.     defs[0].arg1l = -CPI/2.0;
  267.     defs[2].arg1w = MAXLOG*FACT;        /* sqrt */
  268.     defs[2].arg1l = -MAXLOG/2.0*FACT;
  269.     defs[3].arg1w = 2.0*MAXLOG*FACT;    /* exp */
  270.     defs[3].arg1l = -MAXLOG*FACT;
  271.     defs[5].arg1w = 2.0*MAXLOG*FACT;    /* sinh */
  272.     defs[5].arg1l = -MAXLOG*FACT;
  273.     defs[6].arg1w = MAXLOG*FACT;        /* cosh */
  274.     defs[6].arg1l = 0.0;
  275.     defs[8].arg1w = 2.0*MAXLOG*FACT*LOG2E;    /* exp2 */
  276.     defs[8].arg1l = -MAXLOG*FACT*LOG2E;
  277.     defs[9].arg1w = 2.0*MAXLOG*FACT*0.43429448L;    /* exp10 */
  278.     defs[9].arg1l = -MAXLOG*FACT*0.43429448L;
  279.  
  280.  
  281. /* Outer loop, on the test number: */
  282.  
  283.     for (itst = STRTST; defs[itst].name; itst++)
  284.     {
  285.     d = &defs[itst];
  286.     m = 0;
  287.     max = 0.0;
  288.     rmsa = 0.0;
  289.     ave = 0.0;
  290.     fun = d->name;
  291.     ifun = d->inv;
  292.  
  293. /* Absolute error criterion starts with gamma function
  294.  * (put all such at end of table)
  295.  */
  296.     if (d->tstyp == GAMMA)
  297.         printf("Absolute error criterion (but relative if >1):\n");
  298.  
  299. /* Smaller number of trials for Wronksians
  300.  * (put them at end of list)
  301.  */
  302.     if (d->tstyp == WRONK1)
  303.     {
  304.         ntr = WTRIALS;
  305.         printf("Absolute error and only %d trials:\n", ntr);
  306.     }
  307.  
  308.     printf(headrs[d->tstyp], d->nam2, d->nam1);
  309.     printf("random arguments x in " GFORMAT " " GFORMAT " \n",
  310.         d->arg1l, d->arg1l + d->arg1w);
  311.     if (d->nargs == 2)
  312.         printf("random arguments a in " GFORMAT " " GFORMAT " \n",
  313.         d->arg2l, d->arg2l + d->arg2w);
  314.  
  315.     for (i = 0; i < ntr; i++)
  316.     {
  317.         m++;
  318.  
  319. /* make random number(s) in desired range(s) */
  320.         switch (d->nargs)
  321.         {
  322.  
  323.         default:
  324.             goto illegn;
  325.  
  326.         case 2:
  327.             drand(&a);
  328.             a = d->arg2w * (a - 1.0) + d->arg2l;
  329.             if (d->arg2f & EXPSCAL)
  330.             {
  331.             a = exp(a);
  332.             drand(&y2);
  333.             a -= 1.0e-13 * a * y2;
  334.             }
  335.             if (d->arg2f & INT)
  336.             {
  337.             k = a + 0.25;
  338.             a = k;
  339.             }
  340.  
  341.         case 1:
  342.             drand(&x);
  343.             x = d->arg1w * (x - 1.0) + d->arg1l;
  344.             if (d->arg1f & EXPSCAL)
  345.             {
  346.             x = exp(x);
  347.             drand(&a);
  348.             x += 1.0e-13 * x * a;
  349.             }
  350.         }
  351.  
  352.  
  353. /* compute function under test */
  354.         switch (d->nargs)
  355.         {
  356.         case 1:
  357.             switch (d->tstyp)
  358.             {
  359.             case ELLIP:
  360.                 cy1 = (*(fun)) (x);
  361.                 y2 = (*(fun)) (1.0 - x);
  362.                 y3 = (*(ifun)) (x);
  363.                 y4 = (*(ifun)) (1.0 - x);
  364.                 break;
  365.  
  366.             case GAMMA:
  367.                 break;
  368.  
  369.             default:
  370.                 z = (*(fun)) (x);
  371.                 y = (*(ifun)) (z);
  372.             }
  373.             break;
  374.  
  375.         case 2:
  376.             if (d->arg2f & INT)
  377.             {
  378.             switch (d->tstyp)
  379.             {
  380.                 case WRONK1:
  381.                 cy1 = (*fun) (k, x);    /* jn */
  382.                 y2 = (*fun) (k + 1, x);
  383.                 y3 = (*ifun) (k, x);    /* yn */
  384.                 y4 = (*ifun) (k + 1, x);
  385.                 break;
  386.  
  387.                 case WRONK2:
  388.                 cy1 = (*fun) (a, x);    /* iv */
  389.                 y2 = (*fun) (a + 1.0, x);
  390.                 y3 = (*ifun) (k, x);    /* kn */
  391.                 y4 = (*ifun) (k + 1, x);
  392.                 break;
  393.  
  394.                 default:
  395.                 z = (*fun) (k, x);
  396.                 y = (*ifun) (k, z);
  397.             }
  398.             }
  399.             else
  400.             {
  401.             if (d->tstyp == POWER)
  402.             {
  403.                 z = (*fun) (x, a);
  404.                 y = (*ifun) (z, 1.0 / a);
  405.             }
  406.             else
  407.             {
  408.                 z = (*fun) (a, x);
  409.                 y = (*ifun) (a, z);
  410.             }
  411.             }
  412.             break;
  413.  
  414.  
  415.         default:
  416.           illegn:
  417.             printf("Illegal nargs= %d", d->nargs);
  418.             exit(1);
  419.         }
  420.  
  421.         switch (d->tstyp)
  422.         {
  423.         case WRONK1:
  424.             e = (y2 * y3 - cy1 * y4) - 2.0 / (CPI * x);    /* Jn, Yn */
  425.             break;
  426.  
  427.         case WRONK2:
  428.             e = (y2 * y3 + cy1 * y4) - 1.0 / x;    /* In, Kn */
  429.             break;
  430.  
  431.         case ELLIP:
  432.             e = (cy1 - y3) * y4 + y3 * y2 - PIO2;
  433.             break;
  434.  
  435.         default:
  436.             e = y - x;
  437.             break;
  438.         }
  439.  
  440.         if (d->ctrl & RELERR)
  441.         e /= x;
  442.         else
  443.         {
  444.         if (fabs(x) > 1.0)
  445.             e /= x;
  446.         }
  447.  
  448.         ave += e;
  449. /* absolute value of error */
  450.         if (e < 0)
  451.         e = -e;
  452.  
  453. /* peak detect the error */
  454.         if (e > max)
  455.         {
  456.         max = e;
  457.  
  458.         if (e > 1.0e-10)
  459.         {
  460.             /* printf("x %.6E z %.6E y %.6E max %.4E\n",
  461.                x, z, y, max); */
  462.             printf("x " P18EFORMAT " z "  P18EFORMAT
  463.                " y " P18EFORMAT " max " P4EFORMAT " \n",
  464.                x, z, y, max);
  465.             if (d->tstyp >= WRONK1)
  466.             {
  467.             printf("cy1 " P4EFORMAT
  468.                 " y2 " P4EFORMAT
  469.                 " y3 " P4EFORMAT
  470.                 " y4 " P4EFORMAT
  471.                 " k %d x " P4EFORMAT "\n",
  472.                    cy1, y2, y3, y4, k, x);
  473.             }
  474.         }
  475.  
  476. /*
  477.     printf("%.8E %.8E %.4E %6ld \n", x, y, max, n);
  478.     printf("%d %.8E %.8E %.4E %6ld \n", k, x, y, max, n);
  479.     printf("%.6E %.6E %.6E %.4E %6ld \n", a, x, y, max, n);
  480.     printf("%.6E %.6E %.6E %.6E %.4E %6ld \n", a, b, x, y, max, n);
  481.     printf("%.4E %.4E %.4E %.4E %.4E %.4E %6ld \n",
  482.         a, b, c, x, y, max, n);
  483. */
  484.         }
  485.  
  486. /* accumulate rms error    */
  487.         e *= 1.0e16;    /* adjust range */
  488.         rmsa += e * e;    /* accumulate the square of the error */
  489.     }
  490.  
  491. /* report after NTRIALS trials */
  492.     rms = 1.0e-16 * sqrt(rmsa / m);
  493.     printf(" max = " P2EFORMAT " rms = " P2EFORMAT "\n", max, rms);
  494.     }                /* loop on itst */
  495.     return(0);
  496. }
  497.